Projekt 3

W pliku banknotes.csv znadjują się dane opisujące obrazy banknotów. Dane powstały poprzez transformatę falową (wavelett transform) zastosowaną do obrazów w skali szarości rozmiaru 400x400 pikseli. Po zastosowaniu transformaty wyliczono cztery charakterystyki liczbowe obrazu - wariancję, skośność, kurtozę oraz entropię.

Za pomocą modelu regresji logistycznej sprawdź czy za pomocą tej metody jesteśmy w stanie dobrze odróżnić banknoty prawdziwe od fałszywych.

  • Zbuduj i zinterpretuj model regresji logistycznej (w razie otrzymania ostrzeżenia od software’u stosownie należy je skomentować, ale się nim nie przejmować).
  • Zbadaj i zinterpretuj jego charakterystyki liczbowe za pomocą macierzy pomyłek. *Wyrysuj krzywą ROC otrzymanego modelu i podaj AUC. Zinterpretuj otrzymane wyniki.
library(tidyverse)
## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'readr' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'forcats' został zbudowany w wersji R 4.2.3
## Warning: pakiet 'lubridate' został zbudowany w wersji R 4.2.3
library(caret)
## Warning: pakiet 'caret' został zbudowany w wersji R 4.2.3
library(ggfortify)
## Warning: pakiet 'ggfortify' został zbudowany w wersji R 4.2.3
library(knitr)
library(kableExtra)
## Warning: pakiet 'kableExtra' został zbudowany w wersji R 4.2.3
library(plotly)
## Warning: pakiet 'plotly' został zbudowany w wersji R 4.2.3
bank <- readr::read_csv('C:/Users/PC/Downloads/banknote.csv', col_names = FALSE)

Przedstawienie danych

Dane zapisane są w następującej postaci:

  • X1 - wariancja,
  • X2 - skośność,
  • X3 - kurtoza,
  • X4 - entropia,
  • X5 - 1 prawdziwe, 0 fałszywe.
bank %>% count(X5)

Widzimy, że zbiór danych banknote nie jest zbilansowany. Posiada więcej fałszywych banknotów niż prawdziwych.

Normalność

Przy użyciu histogramów oraz testu Shapiro-Wilka sprawdzamy rozkłady naszych zmiennych.

wykres <- plot_ly(alpha = 0.7)
wykres <- wykres %>% add_histogram(x = bank$X1, name = "X1")
wykres <- wykres %>% add_histogram(x = bank$X2, name = "X2")
wykres <- wykres %>% add_histogram(x = bank$X3, name = "X3")
wykres <- wykres %>% add_histogram(x = bank$X4, name = "X4")
wykres <- wykres %>% layout(barmode = "stack")

wykres
shapiro.test(bank$X1)
## 
##  Shapiro-Wilk normality test
## 
## data:  bank$X1
## W = 0.982, p-value = 4.686e-12
shapiro.test(bank$X2)
## 
##  Shapiro-Wilk normality test
## 
## data:  bank$X2
## W = 0.97463, p-value = 8.224e-15
shapiro.test(bank$X3)
## 
##  Shapiro-Wilk normality test
## 
## data:  bank$X3
## W = 0.92709, p-value < 2.2e-16
shapiro.test(bank$X4)
## 
##  Shapiro-Wilk normality test
## 
## data:  bank$X4
## W = 0.91488, p-value < 2.2e-16

Z testu Shapiro-Wilka \(p-value < 0,05\), zatem odrzucamy \(H_0\) i stwierdzamy, że wariancja, skośność, kurtoza oraz entropia nie mają rozkładu normalnego.

Model regresji logistycznej

regresja_logistyczna <- glm(X5~X1+X2+X3+X4,data=bank,family = "binomial")
## Warning: glm.fit: dopasowane prawdopodobieństwa numerycznie okazały się być 0
## lub 1
summary(regresja_logistyczna)
## 
## Call:
## glm(formula = X5 ~ X1 + X2 + X3 + X4, family = "binomial", data = bank)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.70001   0.00000   0.00000   0.00029   2.24614  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.3218     1.5589   4.697 2.64e-06 ***
## X1           -7.8593     1.7383  -4.521 6.15e-06 ***
## X2           -4.1910     0.9041  -4.635 3.56e-06 ***
## X3           -5.2874     1.1612  -4.553 5.28e-06 ***
## X4           -0.6053     0.3307  -1.830   0.0672 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1885.122  on 1371  degrees of freedom
## Residual deviance:   49.891  on 1367  degrees of freedom
## AIC: 59.891
## 
## Number of Fisher Scoring iterations: 12
  • Estimate: szacunki przechwycenia \((\beta_0)\) i współczynnika beta związanego z każdą zmienną przewidującą
  • Std.Error: błąd standardowy oszacowań współczynników. Przedstawia on dokładność współczynników. Im większy błąd standardowy, tym mniejsza pewność co do oszacowania.
  • wartość z: statystyka Z, która jest estymacją współczynnika (kolumna 2) podzieloną przez błąd standardowy estymacji (kolumna 3)
  • Pr(>|z|): Wartość p odpowiadająca statystyce z. Im mniejsza wartość p, tym bardziej znaczący jest szacunek.

Różnica pomiędzy odchyleniem zerowym a odchyleniem resztowym mówi nam, że model jest dobrze dopasowany. Większa różnica jest lepsza dla modelu. Null deviance to wartość, gdy w równaniu mamy tylko przechwyt bez żadnych zmiennych, a Residual deviance to wartość, gdy bierzemy pod uwagę wszystkie zmienne. Ma to sens, aby uznać model za dobry, jeśli ta różnica jest wystarczająco duża.

W naszym przypadku wartość funkcji logitowej \(ln\bigg(\frac{p}{1-p}\bigg)\) z autentyczności banknotów dla wariancji, skośności, kurtozy oraz entropii, wynosi \(7,3218\).

  • Różnica funkcji logitowej autentyczności pomiędzy wartością wariancji wynosi \(-7,8593\), otrzymaliśmy \(p-value < 0,05\), więc odrzucamy \(H_0\). Nie istnieją żadne różnice pomiędzy autentycznością, a wariancją.
  • Różnica funkcji logitowej autentyczności pomiędzy wartością skośności wynosi \(-4,1910\), otrzymaliśmy \(p-value < 0,05\), więc odrzucamy \(H_0\). Nie istnieją żadne różnice pomiędzy autentycznością, a skośnością.
  • Różnica funkcji logitowej autentyczności pomiędzy wartością kurtozy wynosi \(-5,2874\), otrzymaliśmy \(p-value < 0,05\), więc odrzucamy \(H_0\). Nie istnieją żadne różnice pomiędzy autentycznością, a kurtozą.
  • Różnica funkcji logitowej autentyczności pomiędzy wartością entropii wynosi \(-0,6053\), otrzymaliśmy \(p-value > 0,05\), więc nie mamy podstaw do odrzucenia \(H_0\), zatem istnieją różnice pomiędzy autentyycznością banknotu, a jego entropią. Co oznacza, że wpływa ona na autentyczność banknotu.
dec_bond <- regresja_logistyczna$coefficients
logistic2_slope = -dec_bond[2]/dec_bond[3]
logistic2_intercept = -dec_bond[1]/dec_bond[3]
ggplot(bank, aes(x = X1, y = X2, color = as.factor(X5))) + 
  geom_point() + 
  geom_abline(slope = logistic2_slope, intercept = logistic2_intercept, color = 'blue', linetype = 'dashed') + 
    labs(color = 'X5')

ggplot(bank, aes(x = X1, y = X3, color = as.factor(X5))) + 
  geom_point() + 
  geom_abline(slope = logistic2_slope, intercept = logistic2_intercept, color = 'blue', linetype = 'dashed') + 
    labs(color = 'X5')

ggplot(bank, aes(x = X1, y = X4, color = as.factor(X5))) + 
  geom_point() + 
  geom_abline(slope = logistic2_slope, intercept = logistic2_intercept, color = 'blue', linetype = 'dashed') + 
    labs(color = 'X5')

Wykresy przedstawiają hiperpłaszyznę decyzyjną naszego modelu.

Walidacja jakości modelu

logit_prediction <- predict(regresja_logistyczna, bank, type = 'response')
head(logit_prediction)
##            1            2            3            4            5            6 
## 2.220446e-16 2.220446e-16 2.185822e-10 2.220446e-16 4.579103e-01 2.220446e-16
logistic2_predictions <- ifelse(logit_prediction > 0.5, 1, 0)
head(logistic2_predictions)
## 1 2 3 4 5 6 
## 0 0 0 0 0 0
library(caret)
caret::confusionMatrix(data = as.factor(logistic2_predictions), reference = as.factor(bank$X5))
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 757   6
##          1   5 604
##                                          
##                Accuracy : 0.992          
##                  95% CI : (0.9857, 0.996)
##     No Information Rate : 0.5554         
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0.9838         
##                                          
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.9934         
##             Specificity : 0.9902         
##          Pos Pred Value : 0.9921         
##          Neg Pred Value : 0.9918         
##              Prevalence : 0.5554         
##          Detection Rate : 0.5517         
##    Detection Prevalence : 0.5561         
##       Balanced Accuracy : 0.9918         
##                                          
##        'Positive' Class : 0              
## 
table <- data.frame(confusionMatrix(as.factor(logistic2_predictions), as.factor(bank$X5))$table)

plotTable <- table %>%
  mutate(goodbad = ifelse(table$Prediction == table$Reference, "good", "bad")) %>%
  group_by(Reference) %>%
  mutate(prop = Freq/sum(Freq))

ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
  geom_tile() +
  geom_text(aes(label = Freq), vjust = .5, fontface  = "bold", alpha = 1) +
  scale_fill_manual(values = c(good = "green", bad = "red")) +
  theme_bw() +
  labs(title = "Macierz modelu regresja_logistyczna") +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlim(rev(levels(table$Reference)))

Przy pomocy biblioteki caret stworzyliśmy macierz pomyłek. Jej charakterystyki liczbowe przedstawiają się w następujący sposób:

  • True-Negative - fałszywe próbki, które zostały sklasyfikowane jako fałszywe (757),
  • True-Positive - prawdziwe próbki, które zostały sklasyfikowane jako prawdziwe (604),
  • False-Negative - prawdziwe próbki, które zostały sklasyfikowane jako fałszywe (5),
  • False-Positive - fałszywe próbki, które zostały sklasyfikowane jako prawdziwe (6).

Dokładność to liczba próbek poprawnie sklasyfikowanych spośród wszystkich próbek obecnych w zbiorze testowym:

\(Accuracy = \frac{TP+TN}{TP+FP+TN+FN}\), u nas wynosi \(0,992\).

\(NIR > 0,5\), zatem nasz model nie jest gorszy od rzutu monetą.

Współczynnik Kappa może być użyty do oceny dokładności klasyfikacji. Przyjmuje wartości \([-1,1]\) wartość większa niż 0 wskazuje, że klasyfikacja jest znacznie lepsza niż losowa.

Czułość to liczba próbek rzeczywiście należących do klasy pozytywnej spośród wszystkich próbek, które zostały przewidziane przez model jako należące do klasy pozytywnej.

\(Sensitivity = \frac{TP}{TP+FN}\).

Swoistość to liczba próbek przewidywanych prawidłowo jako należące do klasy negatywnej spośród wszystkich próbek w zbiorze danych, które rzeczywiście należą do klasy negatywnej.

\(Sensitivity = \frac{TN}{TN+FP}\).

Zbiór treningowy i testowy

library(caret)
train_test_split <- createDataPartition(bank$X5, list = FALSE, p=0.75)
bank_train <- bank[train_test_split,]
bank_test <- bank[-train_test_split,]
cat(dim(bank_train),dim(bank_test))
## 1029 5 343 5
table(bank$X5)
## 
##   0   1 
## 762 610

Wszystkie banknoty prawdziwe i fałszywe dostępne.

table(bank_train$X5)
## 
##   0   1 
## 568 461

Banknoty treningowe, tylko \(\frac{3}{4}\) banknotów. (p=0.75)

table(bank_test$X5)
## 
##   0   1 
## 194 149

Te banknoty, które odjeliśmy od wszystkich banknotów by otrzymać banknoty treningowe.

kable(summary(bank_train), escape = F) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
X1 X2 X3 X4 X5
Min. :-7.0421 Min. :-13.773 Min. :-5.2861 Min. :-8.5482 Min. :0.000
1st Qu.:-1.8439 1st Qu.: -2.125 1st Qu.:-1.4801 1st Qu.:-2.1888 1st Qu.:0.000
Median : 0.3798 Median : 2.084 Median : 0.6766 Median :-0.5724 Median :0.000
Mean : 0.4016 Mean : 1.777 Mean : 1.4926 Mean :-1.1548 Mean :0.448
3rd Qu.: 2.8237 3rd Qu.: 6.777 3rd Qu.: 3.4356 3rd Qu.: 0.3948 3rd Qu.:1.000
Max. : 6.5633 Max. : 12.952 Max. :17.9274 Max. : 2.4495 Max. :1.000 
bank_train_model <- glm(X5 ~ X1+X2+X3+X4 , data = bank_train, family='binomial')
## Warning: glm.fit: dopasowane prawdopodobieństwa numerycznie okazały się być 0
## lub 1
summary(bank_train_model)
## 
## Call:
## glm(formula = X5 ~ X1 + X2 + X3 + X4, family = "binomial", data = bank_train)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.84582   0.00000   0.00000   0.00047   2.04983  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   6.8728     1.7249   3.984 6.76e-05 ***
## X1           -7.2392     1.9281  -3.755 0.000174 ***
## X2           -4.0288     1.0298  -3.912 9.15e-05 ***
## X3           -4.9895     1.3060  -3.820 0.000133 ***
## X4           -0.6681     0.3912  -1.708 0.087669 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1415.35  on 1028  degrees of freedom
## Residual deviance:   35.17  on 1024  degrees of freedom
## AIC: 45.17
## 
## Number of Fisher Scoring iterations: 12
logit_prediction_bank_train_model <- predict(bank_train_model, bank_train, type = 'response')
head(logit_prediction_bank_train_model)
##            1            2            3            4            5            6 
## 2.220446e-16 2.220446e-16 1.762616e-09 2.220446e-16 5.715952e-01 2.220446e-16
logistic2_predictions_bank_train_model <- ifelse(logit_prediction_bank_train_model > 0.5, 1, 0)
head(logistic2_predictions_bank_train_model)
## 1 2 3 4 5 6 
## 0 0 0 0 1 0
library(caret)
caret::confusionMatrix(data = as.factor(logistic2_predictions_bank_train_model), reference = as.factor(bank_train$X5))
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 562   5
##          1   6 456
##                                          
##                Accuracy : 0.9893         
##                  95% CI : (0.981, 0.9947)
##     No Information Rate : 0.552          
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0.9784         
##                                          
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.9894         
##             Specificity : 0.9892         
##          Pos Pred Value : 0.9912         
##          Neg Pred Value : 0.9870         
##              Prevalence : 0.5520         
##          Detection Rate : 0.5462         
##    Detection Prevalence : 0.5510         
##       Balanced Accuracy : 0.9893         
##                                          
##        'Positive' Class : 0              
## 
table <- data.frame(confusionMatrix(as.factor(logistic2_predictions_bank_train_model), as.factor(bank_train$X5))$table)

plotTable <- table %>%
  mutate(goodbad = ifelse(table$Prediction == table$Reference, "good", "bad")) %>%
  group_by(Reference) %>%
  mutate(prop = Freq/sum(Freq))

ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
  geom_tile() +
  geom_text(aes(label = Freq), vjust = .5, fontface  = "bold", alpha = 1) +
  scale_fill_manual(values = c(good = "green", bad = "red")) +
  theme_bw() +
  labs(title = "Macierz modelu bank_train") +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlim(rev(levels(table$Reference)))

fmpreds <- predict(bank_train_model, bank_train, type = 'response')
fmpreds_classes <- ifelse(fmpreds > 0.5, 1, 0)
baseline_cm <- confusionMatrix(factor(fmpreds_classes), factor(bank_train$X5))
baseline_cm
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 562   5
##          1   6 456
##                                          
##                Accuracy : 0.9893         
##                  95% CI : (0.981, 0.9947)
##     No Information Rate : 0.552          
##     P-Value [Acc > NIR] : <2e-16         
##                                          
##                   Kappa : 0.9784         
##                                          
##  Mcnemar's Test P-Value : 1              
##                                          
##             Sensitivity : 0.9894         
##             Specificity : 0.9892         
##          Pos Pred Value : 0.9912         
##          Neg Pred Value : 0.9870         
##              Prevalence : 0.5520         
##          Detection Rate : 0.5462         
##    Detection Prevalence : 0.5510         
##       Balanced Accuracy : 0.9893         
##                                          
##        'Positive' Class : 0              
##

Krzywa ROC

Analiza krzywej charakterystyki operacyjnej odbiornika (ROC — Receiver operating characteristic) jest użytecznym narzędziem oceny dokładności predykcji modelu poprzez wykreślenie czułości wobec (1-swoistości) testu klasyfikacyjnego (przy progu zmieniającym się w całym zakresie wyników testu diagnostycznego).

roc_rl <- pROC::roc(response = bank_train$X5, predictor = fmpreds)
pROC::ggroc(roc_rl, legacy.axes = TRUE) + geom_abline(slope = 1, intercept = 0) +
  labs(title = "Krzywa ROC bank_train")+
  theme(plot.title = element_text(hjust = 0.5))

roc_rl$auc
## Area under the curve: 0.9999

Krzywa kształtem bardzo przypomina trójkąc co oznacza, że nasz model jest prawie “idealny”, a pole pod wykresem wynosi \(0,9998\).

Podsumowanie

Dzięki tej metodzie jesteśmy w stanie odróżnić prawdziwe banknoty od fałszywych, ale zdarza się, że nasz model popełni błąd. Tworząc model treningowy z \(\frac{3}{4}\) wszystkich banknotów uzyskaliśmy lepsze współczynniki No Information Rate, Sensitivity oraz Specificity. Co za tym idzie nasz model jest bliższy modelowi o idealnych parametrach.